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The Galton board is often used in the classroom as a demonstration experiment for the probability distribu- 
tion of independent events, The falling particles show normal distribution in the lower slots. However, when 
many particles fall simultaneously, the influence between the particles cannot be ignored, and the distribution 
function in demonstration experiments usually deviates from the Gaussian distribution function. In this paper, 
we use simulation experiments to study the impact of interactions on distribution laws under different rate of 
flow. We use a simplified model to parameterize a the strength of the interaction, which describes the impact of 
the particle numbers on direction of particles. The simulation experiment found that the flow rate Nsm (number 
of particles falling simultaneously) and the interaction factor a jointly changed the distribution function. After 
introducing the interaction, the Gaussian distribution and the broadened Gaussian distribution, as well as the 
combination of multiple Gaussian distributions, cannot fit the experimental data well. The Fermi Dirac distribu- 
tion unexpectedly conforms to the experiment, which reflects the ’quantum’ effect of the classic toy experiment. 
Particles appear as a whole in a certain place, and the probability is affected by interactions. Probability inter- 
pretation and wholeness are the essence of quantum physics. The exclusivity between particles leads to negative 
Chemical potential in the distribution function, and the temperature parameter increases with the interaction in- 
tensity a and flow rate Nsm. The relations between parameters can be expressed as a conjecture formula within 
large parameters range. 
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I. INTRODUCTION Galton board 
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The Galton board [1] is usually used as a toy experiment to 
demonstrate probability distributions of independent events. 
The simple model also can study metal electron gas systems 
after considering dynamic details [2, 5, 8]. Besides the studies 
of physics, Galton boards are also used in medical and traffic 
research [9, 10]. We use the following simplified model to 
discuss the working principle of Galton board. As shown in 
Fig. 2, the Galton board device consists of evenly spaced ver- 
tical slots below and regularly arranged horizontal bars in the 
middle. The upper container contains a large amount of par- 
ticles. Open the small outlet at the bottom of the container, 
and a large number of particles will flow down from top to 
bottom, colliding with the horizontal bars of each layer and FIG. 1: The imagine of particles randomly falling in a Galton board 
falling into the slots. The Galton board simulates the track of generated by MATLAB. The Galton board consists of 9 slots below 
particle motion, counts the number of particles falling into the and 9-layer horizontal bars(the dots in the figure represent the spaces 


slot and their proportion to the total number of particles. between bars ). Particles fall one by one, The 685th particle enters 
: : 5th slot through a random path as shown in the figure. 
The Galton board has been easily approached using nu- 


merical simulations experiments[4], with the development of 
numerical calculation software. Numerical calculations can 
even calculate the dynamic track of particle falling [6, 7] . cle falls into one of m slots and the distribution of position 
In our simulation experiments, particles move left or right by satisfies the binomial distribution: 

one grid according to the set probability for each layer they 

fall, which is equivalent to one binomial selecting distribution pi = Ci O m! 
event. If particles fall alone or ignore interactions between m gm (i — m)li!2m’ 
particles, the probability of left and right shifts is the same. 
In this case, after passing m-layer horizontal bars, one parti- 


(1) 


Here, PŻ, is the probability of particles falling into the i-th 
slot in total m slots. 


When the number of layers/slots in the Galton board is large 
*Electronic address: mikehh2021 @163.com enough(m > 9), the distribution is approximately a Gaussian 
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FIG. 2: The simulation experiment values and fitted Gaussian distri- 
bution function(standard deviation 0 = 1.45) compare to the bino- 
mial distribution P, as total slots number m = 9. 
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Here, x = (i — m)/2 is the slot position where particles 
entered. The standard deviation 0 = 1.76 when m = 13 
and ø = 1.45 when m = 9 which according to the following 
theoretical formula 
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We simulated this process by MATLAB and obtained re- 
sults that are consistent with theory. 


(3) 


II. INTRODUCTION AND PARAMETERIZATION OF 
INTERPARTICLE SQUEEZING EFFECTS 


Introduction and parameterization of interparticle squeez- 
ing effects. 

The probability of particles moving left or right is influ- 
enced by the number of particles on both sides, with fewer 
particles and more vacancies corresponding to a higher proba- 
bility. To simplify the problem, we assume that the probability 
of selecting left or right moving is influenced by the relative 
difference of particle numbers at two sides. If the number of 
particles in same grid is Nm, on the left is Nz, and Nr on 
the right, then the probability of left shifting will increased by 
Apr = $(Nr—-Nz)/(Ni + Nm + Np), where the factor a 
represents the strength of the interaction, which is determined 
by the ratio of space between horizontal bars to particle diam- 
eter. The total probability of left and right shifting is normal- 
ized as pr, + pr = 1. So, the probabilitis 
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FIG. 3: The imagine of particles randomly falling in a Galton board 
generated by MATLAB. The Galton board consists of 13 slots below 
and 13-layer horizontal bars. Nsm = 5 particles fall simultaneously, 
interaction factor œ = 0.5. The five particle paths are shown in the 
figure. 
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FIG. 4: The fitted Gaussian distribution function(standard deviation 
o = 2.03) compare to the simulation experiment values as total slots 
number m = 17. The consistency of data is better than that in Figure 
2 because the layer number m is larger. when a is zero, the distri- 
bution functions of Nsm = 1 (Left) and Nsm = 10(Right) are the 
same Gaussian function. 


We simulated the process of Nsm particles falling simulta- 
neously by MATLAB as shown in Fig. 3. 

If a is zero, this means the particles are tiny and the inter- 
action between particles is negligible. In this case, no mat- 
ter how many particles fall simultaneously, the probability 
pL = pr = 1/2 and the distribution function is the same 
as that of free particle passing through the Galton board. The 
simulation experiment results verify this conclusion, as shown 
in Fig. 4, the simultaneous falling particle number Nsm = 1 
and Nsm = 5, the distribution functions are the same Gaus- 
sian function. 

When both Nsm and a increase, the simulation distribution 
obviously deviates from the Gaussian function and can not 
fitted by a broaden Gaussian function as shown in Fig. 5. 

Firstly, it is assumed that the result of the simulation exper- 
iment is two extended Gaussian distributions with left/right 
shift in the mean value. We use two symmetric Gaussian func- 
tions with the center value shift u to fit the simulation experi- 
ment results 
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FIG. 5: The simulation experiment values and best-fitting Gaussian 
distribution function(standard deviation o = 2.7), Nsm = 10 and 
a = 0.5. The simulation distribution obviously deviates from the 
Gaussian function. 
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FIG. 6: The simulation experiment values and best-fitting shifted 
Gaussian functions(standard deviation 0 = 2.25, mean value pp = 
1.56), Nsm = 20 and a = 0.5 (Left). (@ = 2.27, mean value 
H = 1.96), for Nsm = 7 and a = 1(Right) 


Then, two shift Gaussian functions has a much better fit with 
than simulation results as shown in Fig. 6 when Nsm is large 
and a not close to 1. But two Gaussian fitting is not good 
when small Nsm and a > 0.7. 

Meanwhile, the best-fitting parameters u and o are pre- 
sented as complex function curves with the different Nsm and 
a as shown in Fig. 7. 

The u and o values of the Gaussian function increase with 
Nsm and a parameters. As shown in the Fig. 4. The re- 
sults indicates that the number of particles falling simultane- 
ously Nsm and the interaction factor œ cause the particles to 
disperse on both sides due to collision repulsion during the 
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FIG. 7: The best-fitting parameters js and o vs variant a for different 
Nsm 
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FIG. 8: The simulation experiment with fixed a = 0.5 and variant 
Nsm = 5,10, 1000, compared with the best fitting Fermi Dirac dis- 
tributions 


falling process. 
In summary, the extended Gaussian distribution cannot de- 
scribe the experimental results. 


A. Simulation Models 


Considering the exclusivity caused by interactions and 
wholeness of particles, we attempt to fit the experiment us- 
ing the Fermi Dirac distribution function: 


1 x? 
ae ae (7) 
e i> +1 


f(x) x 


The temperature parameter 1 + b and Fermi energy(Chemical 
potential) Ær in this function are obtained through fitting ex- 
periments. The exclusivity between particles leads to negative 
Chemical potential when Nsm is larger, and the temperature 
parameter increases with the interaction intensity œ and flow 
rate Nsm. Fermi Dirac distributions functions can fit experi- 
ments well within the full parameter range for Nsm = 2 — oo 
and a = 0 — 1. The best fitting Fermi Dirac distribution 
curves compared with experiments are shown in Fig. 9 and 
Fig. ??. 

Although the Fermi Dirac distribution function can describe 
experiments well, the fitting parameters especially Fermi en- 
ergy Er have significant fluctuations as shown in Fig. 10. The 
reason for the fluctuations is that most particles are distributed 
in the 3Sigma range about 7 sluts, which leads to the avail- 
able data for fitting is not enough. 

In order to find the relationship between the Fermi Dirac 
distribution parameters b, — Ep and experimental parameters 
a and Nsm, the more slut-number m is needed. When the 
slut-number m = 25, the fluctuation of b parameter is re- 
duced significantly. The parameter b exhibits a simple lin- 
ear dependence on the interaction intensity a, when flow rate 
Nsm > 10. 

The ’temperature’ parameter b can be simply expressed as: 


b = kae 7°; (8) 
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FIG. 9: The simulation experiment with variant œa and Nsm, com- 
pared with the best fitting Fermi Dirac distributions 


FIG. 10: The the best fitting parameters b and — Ep with variant a 
and Nsm, Fermi energy Er have significant fluctuations as m = 17. 


Here, the slop ky slowly increases with the flow rate Nsm 
as the conjecture formula 
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(9) 
The exponential factor 7, rapidly decreases to zero with in- 
creasing Nsm 
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When Nsm is large, Fermi energy Ep is not sensitive to 
the quality of fitting. Otherwise, the fitted Fermi energy Ep 


FIG. 11: The the best fitting parameters b and — Ep with variant a 
and Nsm as m = 25, the parameter b can be simply expressed as a 
formula b = Kyae "°° . 


à 


N 


FIG. 12: The parameter € vs variant «œ and different Nsm, the curves 
are simple linear functions as parameter Nsm > 10. When Nsm < 
10, € can be simply expressed as a formula € = Keae"’ . 
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FIG. 13: The simple linear relationship conjecture formulas can well 
fit the experiment when the flow rate Nsm is large. 


is not accurate and has large fluctuations as shown in Fig. 11. 
It’s need to find a new parameter that is sensitive to fitting. We 
find that the ratio of the FD function’s second derivative to the 
value of function at the origin is a good parameter. The new 
sensitive parameter € (the ratio minus 1 )is defined as 


fab kia (11) 


The curves of € are simple linear functions as parameter 
Nsm > 10 as shown in Fig. 12. When Nsm < 10, € can be 
simply expressed as a formula 


E = Kea" (12) 


When the flow rate Nsm is large, these curves are simple 
linear(7j7~¢ = 0) to a and the slop ke ~ 2.5. Therefore, pa- 
rameters b and € can be determined using above simple linear 
functions of œ for Nsm > 10. The values of Ep can be ob- 
tained using conjectured b and € as 
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(13) 

The above simple linear relationship conjecture formulas 
can well fit the experiment when the flow rate Nsm > 10 as 
shown in Fig. 13. 

There are still some fluctuations in the curve when Nsm < 
10 , which come from uncertain fitting Æp. In order to de- 
termine the slope K¢ and exponential factor ne when Nsm is 
small , we need to eliminate these fluctuations as much as pos- 
sible. So, we fixed the b parameters using the above conjecture 
formula, fitted the experiment data again to obtain the smooth 
curves of Ep and £. 


FIG. 14: The parameter € vs variant a and different small Nsm, the 
€ curves can be simply expressed as a formula € = keae"é* . The 
best-fitting Er also become smooth and conform to the conjecture 
formula. 
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FIG. 15: The non-linear conjecture formulas can well fit the experi- 
ment when the flow rate Nsm is small. 


It can be seen in Fig. 14 that the € curves can be simply 
expressed as a formula € = keae”** . The slop srg increases 
with the flow rate Nsm from 2 to 10 as the conjecture formula 
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The exponential factor 7¢ is small and changes according 
to the following polynomial of 1/Nsm 


ne = 7.687 — 14.58? (15) 


In the Fig. 15, it can be seen that the above conjecture for- 
mula for small Nsm can well fit the experiment data. 

Now, we can use a similar conjecture formula to express the 
distribution of m = 17. The parameter b also expressed as: 


b = Kypae 7%; (16) 


The exponential factor 7, rapidly decreases to zero with in- 
creasing Nsm 
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When Nsm > m/2, the parameter b is still linear to œ and 
the slop &» slowly increases with the flow rate Nsm as the 
conjecture formula 


m— 1 
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(18) 


There are only one coefficient changed from 1.99 for m = 
25 to 1.68 for m = 17. 
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FIG. 16: The simple linear relationship conjecture formulas € = 
2.6a and b = 1.68(1 — *)a can well fit the experiment when the 
flow rate Nsm > m/2. 
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FIG. 17: The parameter £ vs variant a and different small Nsm, the 
€ curves can be simply expressed as a formula € = Keae”é“ . The 
best-fitting Ær also become smooth and conform to the conjecture 
formula. 


When Nsm > m/2, the parameter £ is still linear to œ and 
the slop «g is a constant 2.6 which is a little fine-tuning. 

We fixed the b parameters using the above conjecture for- 
mula, fitted the experiment data for m = 17 again to obtain 
the smooth curves of Ep and £. 

When Nsm < m/2, the kp become tiny value and the cur- 
vature of € must be considered as shown in Fig. 17. € still can 
be expressed as € = gae”: ~, Here, kg becomes 


1 
Ke = (2.17 + 2.368e! 532) kp; B= N-I (19) 


The exponential factor 7¢ changes according to the follow- 
ing polynomial of 1/Nsm 


ne = 128° — 226°. (20) 


In Fig. 18, it can be seen that the above conjecture formulas 
for small Nsm can well fit the experiment data m = 17. 
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FIG. 18: The non-linear conjecture formulas for m = 17 can well fit 
the experiment when the flow rate Nsm is small. 
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FIG. 19: The figure demonstrates the dependence of —EF /(1 + b) 
on experimental parameters a and Nsm for m = 17. The Fermi 
Dirac distribution will degenerate into a Gaussian distribution func- 
tion when a is small. The upper limit on a for degradation will 
increases with Nsm. 
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FIG. 20: When the value of m = 25, — Ep /(1 + b) will increase. 
The upper limit on a for degradation even becomes 0.5 when Nsm 
is large. 


The Fermi Dirac distribution will degenerate into a Gaus- 
sian distribution function when Ep > 2.5. 

The Fig. 19 demonstrates the dependence of —E/(1 + b) 
on experimental parameters a and Nsm form = 17. The 
Fermi Dirac distribution will degenerate into a Gaussian dis- 
tribution function when œ is small. The upper limit on a 
for degradation will increases with Nsm. When Nsm reaches 
1000, The distribution will degenerate as a < 0.3. 

The increase of the sluts number m leads to the growth of 
—Ep/(1 + b), making degradation more likely to occur. It 
can been seen from Fig. 20 that the distribution even degen- 
erates when a < 0.5 if m increased to 25. This phenomenon 
indicates that m represents the number of particle states, and 
when the number of states is large, the states is almost con- 
tinuous spectrum. The quantum distribution will degenerate 
to the Classical limit, and the Fermi Dirac distribution will 
become the Classical Boltzmann distribution which is equiv- 
alent to Gaussian distribution. In order to verify the Classical 
limit of large m, we directly extend m to 101. The conjecture 
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FIG. 21: When m = 101, the value of —Ep/(1 + b) is large in 
most range. Regardless of the magnitude of the interaction, the dis- 
tribution function can degenerate into a Gaussian distribution when 
Nsm > 100. 


formula is still valid, just adjust the «x, coefficient. 
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; (21) 
while m = 101. 

Meanwhile, when Nsm > 20, the parameter £ is still linear 
to a and the slop «g is a constant 2.5 which is same with 
m = 25; 

The large values of —E-/(1 + b) in Fig. 21 indicate that 
classical correspondence does occur in the case of continuous 
spectrum. The distribution function returns to a broadened 
Gaussian distribution. And the broadening ratio of sigma? is 
a simple linear relationship with interaction factor a. 


II. REAL GALTON BOARD EXPERIMENT 


By now, all the study is based on simulation experiment 
data. In order to verify whether the above analysis results 
match the real Galton Board experiment, a real instrument 
with m=17 is used to obtain real data. It is difficult to measure 
the interaction intensity of particles in real Galton Board. The 
flow rate of particles falling simultaneously is also difficult to 
accurately control. So the values of a and N,,,, can only be 
determined based on experimental data. 

The real Galton Board is much more complex than ideal 
model. Particles falling alone do not follow a binomial distri- 
bution, because the step size of particle movement is not fixed 
+0.5. The Gaussian distribution is still regarded as the prob- 
ability distribution function when particles fall alone, but the 
standard deviation o # 2.02 and needs to be measured. We 
first use small particles flow to obtain the standard deviation. 
The true standard deviation is 2.3 measured from the particle 
distribution data in Fig. 22. 

Real experimental devices can only control the flow rate 
within a rough range, so Nsm is not fixed. The particles in- 
teraction factor a also cannot be determined as a constant in- 
dependently of Nsm , just like the ideal model. Therefore, 
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FIG. 22: The above figure shows the distribution of small flow real 
experiments and the following figure shows the high flow. 
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FIG. 23: The Fermi Dirac distribution function with b = 1.15 and 
— Ep = 3 is consistent with the real experimental distribution data. 


we can only roughly determine the interaction factor œ corre- 
sponding to the flow rate Nsm within a certain range from the 
experimental results. 

The basic standard deviation measured for this experimen- 
tal device is o = 2.3, this standard deviation corresponds to 
the value m = 22 in the ideal model.According to the conclu- 
sion of the simulation experiment, the parameter b is linear to 
a, When Nsm > m/2, and the slop kẹ slowly increases with 
the flow rate Nsm as the conjecture formula for m = 22 
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Based on the real experimental results of Fig. 23, the flow 
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FIG. 24: The Bose-Einstein distribution function(Right) fits the ex- 
periment much better than Gaussian distribution(Left) 
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FIG. 25: The distribution width 1 + b decreases with the increase of 
—aand Nsm, indicating that attraction and particle density enhance 
Bose Einstein condensation 


rate Nsm is approximately 20-50. According to the conjecture 
formula of the simulation experiment, the equivalent interac- 
tion factor of this device a ~ 0.6—0.7 when Nsm œ 20—50. 


IV. ATTRACTION AND BOSE-EINSTEIN DISTRIBUTION 


The interesting results of simulation experiments of Gal- 
ton boards with exclusivity between particles can be easily 
extended to other hypothetical models. It is possible to imag- 
ine the situation where particles have attraction interaction. 
We can conveniently set a to negative values to parameterize 
the Toy model. Like the exclusive Galton board, when the 
attraction is in the middle, the probability distribution devi- 
ates from the Gaussian distribution which can been seen in 
Fig. 24(Left). It’s reasonable to assume that particles distribu- 
tion law would be Bose-Einstein distribution functions. From 
Fig. 24(Right), it can be seen that the Bose-Einstein distribu- 
tion function fits the experiment much better. 

We continue to use 1 + b as the temperature parameter in 
the distribution functions, and jz is the Chemical potential. At- 
traction will cause b to be negative, while its value close to -1 
representing the appearance of Bose-Einstein condensation. 


(23) 


The distribution width 1 + b is very narrow when attrac- 
tion —qa is strong and particle number Nsm is large as shown 
in Fig. 25, meaning almost particles condense into the same 
state. 

The distribution square-width 1 + b(also the temperature 
parameter) decreases with —a rapidly and linearly until 1 + 
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FIG. 26: The distribution width 1 + b decreases with —a, and 
—a/(1 + b) is not large. 


b < 0.2 which meaning the particles are only distributed in 
three positions. The parameter Chemical potential — u over 
”temperature” (kT = 1 + b) represents the difference be- 
tween quantum distribution and classical distribution func- 
tions. From Fig. 26(Right), it can be seen that Bose-Einstein 
distributions cannot degenerate into Gaussian distributions 
when —o is in the middle value. 


V. DISCUSSION AND CONCLUSIONS 


Through simulation and real experiments, it was found that 
the distribution of particles in Galton Boards cannot be sim- 
ply represented by a Gaussian distribution or even a binomial 
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